module modStresses

    contains
    
    
    subroutine convertStressVector2Matrix(StressVec, StressMat)
        implicit none
        REAL(8), INTENT(IN) :: StressVec(6)
        REAL(8), INTENT(OUT) :: StressMat(3,3)
        
        ! Local Variables
        INTEGER(4) :: ii
        
        ! Normal Stress {Sigma11, Sigma22, Sigma33
        Do ii = 1, 3    
            StressMat(ii,ii) = StressVec(ii)       
        enddo
        
        ! Shear Stresses
        StressMat(1,2) = StressVec(4)       ! Sigma12
        StressMat(1,3) = StressVec(5)       ! Sigma13
        StressMat(2,3) = StressVec(6)       ! Sigma23
        ! Transpose
        StressMat(2,1) = StressMat(1,2)
        StressMat(3,1) = StressMat(1,3)
        StressMat(3,2) = StressMat(2,3)
    
        return
    end subroutine convertStressVector2Matrix
    
    
    ! Calculates the First Invariant of stress, using the stress vector
    ! Input:
    ! stressVector = s1, s2, s3, s4, s5, s6  where first 3 is Normal stress, and last 3 is Shear Stress.
    ! Output:
    ! firstInvariant = Sum of three normal stress = s1 + s2 + s3
    SUBROUTINE calcFirstInvariant(stressVector, firstInvariant)
        implicit none
        REAL(8), INTENT(IN) :: stressVector(6)
        REAL(8), INTENT(OUT) :: firstInvariant
    
        firstInvariant = stressVector(1) + stressVector(2) + stressVector(3)
        ! or it can also be written
        firstInvariant = SUM(stressVector(1:3))
        
    end subroutine calcFirstInvariant



    ! Calculates the Second Invariant of stress, 
    ! Input:
    ! stressVector = s1, s2, s3, s4, s5, s6  where first 3 is Normal stress, and last 3 is Shear Stress.
    ! firstInvariant = Sum of three normal stress = s1 + s2 + s3
    ! Output:
    ! deviatoricStress - matrix of [3,3], see Eqn 2.48
    ! secondInvariant - scalar, calculated using Eqn 2.51b
    subroutine calcSecondInvariant(stressMatrix, firstInvariant, deviatoricStress, secondInvariant)
        implicit none
        REAL(8), INTENT(IN) :: stressMatrix(3,3)
        REAL(8), INTENT(IN) :: firstInvariant
        REAL(8), INTENT(OUT) :: deviatoricStress(3,3)
        REAL(8), INTENT(OUT) :: secondInvariant
        
        ! Local Variables
        REAL(8) :: UnitMatrix(3,3)
        INTEGER(4) :: ii
        
        ! Create Unit Matrix
        UnitMatrix(:,:) = 0.0d0             ! Initialize Entire Matrix
        Do ii = 1, 3
            UnitMatrix(ii,ii) = 1.0d0       ! Create diagonal of Unit Matrix
        enddo
        
        ! Calculate Deviatoric Component of Stress Tensor - See Eqn 2.48
        deviatoricStress(:,:) = stressMatrix(:,:) - (firstInvariant/3.0d0) * UnitMatrix(:,:)
        
        ! Calculate the second invariant
        secondInvariant = SUM( stressMatrix * stressMatrix )    ! let X = s11*s11 + s12*s12 + s13*s13 + s21*s21 + ....
        secondInvariant = SQRT(0.5d0 * secondInvariant)         ! SQRT(X) / SQRT(2)
        
        RETURN        
    end subroutine calcSecondInvariant       


    subroutine calcThirdInvariant()
        implicit none
        
        ! COMING SOON
        
    end subroutine calcThirdInvariant       


end module modStresses